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Abstract. Ultrasonic guided waves have proved useful for structural health monitoring (SHM) and nondestructive 
evaluation (NDE) due to their ability to propagate long distances with less energy loss compared to bulk waves and due 
to their sensitivity to small defects in the structure. Analysis of actively transmitted ultrasonic signals has long been used 
to detect and assess damage. However, there remain many challenging tasks for guided wave based SHM due to the 
complexity involved with propagating guided waves, especially in the case of composite materials. The multimodal 
nature of the ultrasonic guided waves complicates the related damage analysis. This paper presents results from parallel 
3D elastodynamic finite integration technique (EFIT) simulations used to acquire 3D wave motion in the subject 
laminated carbon fiber reinforced polymer composites. The acquired 3D wave motion is then analyzed by frequency- 
wavenumber analysis to study the wave propagation and interaction in the composite laminate. The frequency- 
wavenumber analysis enables the study of individual modes and visualization of mode conversion. Delamination damage 
has been incorporated into the EFIT model to generate “damaged” data. The potential for damage detection in laminated 
composites is discussed in the end. 
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INTRODUCTION 

Current composite “safety margin” design comes with a weight penalty which would be mitigated by the 
development of accurate and efficient structural health monitoring (SHM) methods to locate and quantify the unique 
damage types that can occur in composite materials (e.g., microcracking and delaminations). Unexpected damage 
can occur in aerospace composites due to impact events or due to stressing of the material during off-nominal 
loading events. The problem of how best to detect such material changes in composite components is still being 
solved, and a solution to this problem will be imperative for safe and functional optimally-designed next-generation 
composite aircraft. 

Ultrasonic guided waves have proved useful for structural health monitoring (SHM) and nondestructive 
evaluation (NDE) due to their ability to propagate long distances with less energy loss compared to bulk waves and 
due to their sensitivity to small defects in the structure [1-5]. Analysis of actively transmitted ultrasonic signals is a 
conventional NDE methodology that has long been used to detect and assess damage in aerospace components. 
Many researchers have contributed to the study of SHM using guided wave propagation methods [2-10]. Although 
recent advances in guided wave based SHM technology have demonstrated the feasibility of detecting and locating 
damage in composite structural components, such as [11-13], there remain many challenging tasks due to the 
complexity involved with propagating guided waves. 

In this paper we present our studies on guided wave propagation in anisotropic carbon fiber reinforced polymer 
(CFRP) materials using wavenumber based data analysis methods. In our study, the waveforms are acquired from 
three-dimensional (3D) elastodynamic finite integration technique (EFIT) that generates both in-plane and out-of- 
plane velocity components. Utilizing wavenumber processing methods, wave interaction with delamination damage 
is studied. The two-dimensional Fourier transform (2D FT) is applied to time-space wavefield data such that they 
are converted to the frequency-wavenumber domain where the guided wave propagation characteristics can be 
assessed more intuitively. A short space 2D FT is also developed and applied to generate a space-wavenumber 
relationship. The wave interaction and mode conversion at the delamination region are revealed using these 
wavenumber processing methods. 

FREQUENCY-WAVENUMBER ANALYSIS 

The Fourier transform has proved to be an enormously useful tool for analyzing time series. The concept of 
Fourier analysis is straightforwardly extended to multidimensional signals such as the time -space wavefield, which 


is in terms of the time variable t and the space variable x. Frequency-wavenumber representation can be obtained by 
applying the 2D FT to the time-space wave field u(t, x ), as [14-16] 

U (/, k) - F 2D [u(t, x)] = | | u(t , x)e~ J ^ ft ~ kx " > dtdx ( 1 ) 

where U(f. ] k) is the frequency-wavenumber spectrum which can be interpreted as an alternative representation of the 
time-space wavefield u(t, x) in terms of the frequency variable / and the wavenumber k. 

Figure 1 gives an example of frequency- wavenumber analysis of a simulated harmonic time-space wavefield u 0 (t, 
x) with a constant unit amplitude (Figure 1(a)). Figure 1(b) gives the frequency- wavenumber spectrum in terms of 
amplitude (color map) versus frequency and wavenumber. The spectrum has a maximum value at (150kHz, 0.68 
rad/mm), indicating the harmonic (frequency, wavenumber) component. Note that the small ripples around the 
spectrum center (150kHz, 0.68 rad/mm) are caused by Fourier transform leakage. 


(a) 

FIGURE 1 2D FT of a single frequency single mode numerical wavefield u(t,x ) : (a) time-space wavefield; (b) 

frequency-wavenumber spectrum. 



SPACE-FREQUENCY-WAVENUMBER ANALYSIS 

The frequency- wavenumber analysis unveils wave propagation characteristics that cannot be explicitly seen in 
the time-space domain. However, this approach does not retain the spatial information. Through a spatial windowing 
technique, spatial information can be retained and the relationship between location and frequency/wavenumber can 
be revealed. 

Similar to the idea of the short time Fourier transform which uses a windowing technique [17], a short space 2D 
FT is developed to perform the space-frequency-wavenumber analysis which can retain the spatial information. By 
this means, frequency- wavenumber spectra at various spatial locations are obtained, resulting in a three-dimensional 
space-frequency- wavenumber spectrum of the wavefield. The short space 2D FT of a wavefield u(t, x) is expressed 
as: 


S(a,f,k ) = P f' u{t,x)W\t,x-a)e K2Kft - kc) dtdx (2) 
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where a is the space index and W(t, x) is the 2D window function. In our study, a Hanning function is used to 
construct the 2D window W(t, x), giving: 
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where D x is the window length in the space dimension. The 2D FT of the windowed portion will yield a frequency- 
wavenumber spectrum, which only contains the local information of the windowed region. By sliding the window 
along the space dimension, the resulting space-frequency- wavenumber spectrum can indicate how the frequency- 
wavenumber components vary in space. 

Figure 2 shows an example of the short space 2D FT for the single frequency single mode time-space wavefield 
u 0 (t, x) in Figure 1(a). An illustration of the moving window is shown in Figure 2(a). By sliding the window along 
the space dimension and stacking the corresponding frequency-wavenumber spectra, a space-frequency- 



wavenumber spectrum is obtained. Figure 2(b) plots the resulting space-wavenumber spectrum at a selected 
frequency of 150 kHz. The spectrum shows wavenumber as a function of propagation distance. 



FIGURE 2 Short space 2D FT of a single frequency single mode numerical wavefield u 0 (t, x ): (a) time-space 
window; (b) space-wavenumber spectrum at the frequency 150kHz. 

EFIT WAVE MOTION SIMULATION 


3D EFIT simulations provided simulated guided wave data for the studies presented in this paper. Previous work 
by the authors has demonstrated how 3D simulated guided wave data can aid in the study of guided wave interaction 
with damage in isotropic and anisotropic materials [16, 18, 19]. Since the anisotropic EFIT equations are quite 
lengthy they are not repeated here, but can be found in prior work by Halkjaer [20]. The EFIT code implemented 
for this paper is custom parallelized code that runs efficiently on cluster computing resources. Additional details 
about the anisotropic EFIT simulation code, including results of code validation studies can be found in prior work 
by one of the authors [19]. 


Wave simulation setup 

Wave propagation in an 8-ply IM7/8552 composite laminate with layup [0 2 /90 2 ] s was simulated for two cases: 1) 
an undamaged composite laminate, 2) a composite laminate with a delamination (void) damage located between 
uppermost 0/90 plies (closest to the plate surface where the incident wave excitation is applied). The specimen 
layout is given in Figure 3. The composite layup was chosen to match composite samples that were fabricated at 
NASA LaRC and will be used in future work. In order to isolate mode conversion behavior from more complex 
scattering and geometric effects, long rectangular delaminations of varying lengths were inserted into the simulation 
for the delamination studies. In future work the authors plan to insert square delaminations to match the shape, size, 
and location of Teflon inserts in a NASA LaRC composite sample. This approach allows for the step-by-step 
establishment of a clear understanding of guided wave behavior and wavenumber analysis techniques for composite 
materials containing delaminations. In future work it will be important to move on to study more realistic damage 
geometries. 

The EFIT code incorporates ply-level material properties, and the lay up is simulated by inputting the 
appropriately rotated ply level material properties for each layer. The material properties (single ply layer) used in 
the EFIT simulations are shown in Table 1, where Ei is in the fiber direction. In order to allow for two spatial steps 
per ply, the spatial step size of the EFIT simulation was set smaller than the minimal requirement of Ax < A min / 1 0 to 
Ax=53.3pm (~'k m J 69 based on dispersion curves for the pristine laminate). For the damaged composite specimen 
the 53.3 pm step size is equivalent to ~A aT1 i n /3 1 in the thinned region above the delamination. The wavelength to step 
size ratio above the delamination is expected to be adequate for accurately calculating wave behavior based on the 
results of prior validation comparisons where a similar ratio was used [19]. The total thickness of the simulated 
laminate is 0.854 mm. Additionally, we note that the temporal step size (4.1 ns) meets stability requirements [19]. 
The delamination was incorporated by implementing stress free boundary conditions at the damage location. The 
simulation incident wave was a 300 kHz 3 -cycle Hanning windowed sine wave. The excitation was introduced to 
the simulation space as a 7 mm diameter ring-type excitation on the top laminate surface in order to approximate the 
behavior of a 7 mm piezoelectric wafer [18]. 


TABLE 1. Material properties of the lamina [21] 


p (kg/m 3 ) 

E t (GPa) 

E 2 (GPa) 

E, (GPa) G,, (GPa) G,, (GPa) v 12 v,, 

1570 

171.4 

9.08 

9.08 5.29 2.80 0.32 0.32 
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FIGURE 3 The specimen layout (unit: mm): (a) top view showing PZT, (b) specimen layup and delamination depth. 




FIGURE 4 Theoretical dispersion curves from DISPERSE™: (a) frequency-wavenumber dispersion curves (90 deg); (b) 
group velocity dispersion curves (90 deg). 




(a) (b) (c) (d) 

FIGURE 5 EFIT simulation out-of-plane (OOP) results at time t=24.3ps on the for the top surface of the composite plate: 
(a) pristine plate; (b) with 5mm delamination; (c) with 10mm delamination; (d) with 20mm delamination. 





Simulation results 


Full wavefield velocity data (in-plane v x and v y components and out-of-plane v z component) were output from the 
EFIT simulation for the entire top surface of the simulated composite sample. The surface velocity data was output 
every 20 time steps in order to reduce the size of the resulting data set. This time output meets sampling 
requirements of the data analysis techniques. 

Dispersion curves for the [0 2 /90 2 ] s were calculated using DISPERSE software , and results for the 90 deg fiber 
direction (corresponding to the point-line in Figure 3 and the results in the following section) are shown in Figure 4. 
At the selected frequency-thickness of the simulation, two guided wave modes are expected to exist. 

Figure 5 shows time domain wavefield images at a single point in time from EFIT simulations of the pristine 
sample and for delamination damage of increasing length, as indicated in the figure caption. The A 0 mode is clearly 
observed from the out-of-plane wavefield images. The faster S 0 mode has much lower out-of-plane amplitude and is 
not readily visible in the images. Wave interactions with the delamination regions can be observed in the three 
damage scenarios, however, specific behaviors such as mode conversion and wavenumber shifts, are not clear from 
the time domain wavefield images. 

WAVE INTERACTION WITH DELAMINATION 


In this section, the 3D in-plane and out-of-plane wave motion data generated by EFIT model are analyzed with 
the frequency-wavenumber and space-frequency-wavenumber approaches in order to understand the wave 
propagation characteristics and wave interaction with the delamination damage. The wave propagation along the 90° 
direction (the point line in Error! Reference source not found.), which passes through the delamination area, is 
studied. Since the v x component is nearly zero, here we show only the in-plane v y and out-of-plane v z components. 
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FIGURE 6 In-plane v y results: (a), (b), (c) and (d) are the time-space wavefields for the pristine plate, plate with 
Z=5mm delamination, plate with Z=10mm delamination and plate with Z=20mm delamination; (e), (f), (g) and (h) are 
the frequency wavenumber spectra corresponding to the time-space wavefields in (a), (b), (c) and (d), respectively; (i), 
(j), (k) and (1) are the space-wavenumber spectra corresponding to the wavefields in (a), (b), (c) and (d), respectively. 
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In-plane v y component 


The time-space wavefields of the in-plane v y component for pristine and delamination cases containing different 
delamination sizes (L= 5, 10 and 20 mm) are shown in Figure 6(a)-(d). The wavefield is plotted as amplitude (in 
color) versus time t and propagating distance x with respect to the actuator position. A fast S 0 mode and a slow A 0 
mode are observed in the pristine plate. The in-plane amplitude of S 0 mode is higher than that of A 0 mode, as 
expected. For the delamination cases, it can be observed that the waves propagate forward and backward within the 
delamination region and the energy is somewhat “trapped” but slowly “leaking” from the damage region. This 
observation is consistent with what has been reported previously in the literature [13]. Although not readily 
apparent with the color scale used in figure 6, we also note that as the S 0 mode impinges and reflects from the 
delamination edge, an A 0 mode is created due to mode conversion. 

In order to leverage frequency-wavenumber information, both 2D FT and short-space 2D FT analyses are applied 
to the data. A window size of 50 mm (twice the wavelength of the S 0 mode in the full plate thickness) was used for 
the short-space 2D FT, in order to obtain good space-wavenumber resolution. Using the 2D FT, the time space 
wavefield is transformed into the frequency-wavenumber domain. The frequency-wavenumber spectra of the 
wavefields for each case are given in Figure 6(e)-(h). The solid black line and red dashed line are theoretical 
dispersion curves. In the spectrum of the pristine plate, both A 0 and S 0 mode can be seen, with S 0 being stronger for 
the in-plane data. The predicted modes are consistent with the time-space wavefield observations. In the 
delamination cases additional wavenumber components (both positive and negative components, representing 
forward and backward wave propagation) appear. In Figure 6(h), we can clearly observe a new frequency- 
wavenumber component which may be related to the waves trapped inside the delamination region. 




FIGURE 7 Out-of-plane v z results: (a), (b), (c) and (d) are the time-space wavefields for the pristine plate, plate with 
L=5mm delamination, plate with Z=10mm delamination and plate with Z=20mm delamination; (e), (f), (g) and (h) are 
the frequency wavenumber spectra corresponding to the time-space wavefields in (a), (b), (c) and (d), respectively; (i), 
(j), (k) and (1) are the space-wavenumber spectra corresponding to the wavefields in (a), (b), (c) and (d), respectively. 





Using short- space 2D FT, the space information can be retained, in addition to the wavenumber information. 
Figures 6(i)-(l) present space-wavenumber spectra at 300 kHz. As shown in Figure 6(1), within the delamination 
region, there exist new wavenumber components (again including both positive and negative wavenumbers 
representing forward and backward wave propagation). The center of the new wavenumber component is around 
0.75 rad/mm which is close to the theoretical wavenumber (0.78 rad/mm) for the S 0 mode in the upmost two ply 
layers (i.e., above the delamination). 


Out-of-plane v z component 

The time-space wavefields of out-of-plane v z component for the pristine and delamination cases (with 
delamination sizes L= 5, 10 and 20 mm) are shown in Figure 7(a)-(d). A strong A 0 mode can be clearly seen in the 
wavefield of pristine plate in 

Figure 7(a). In contrast, the out-of-plane amplitude of S 0 mode is very weak. At the delamination region, wave 
interactions with the delamination appear. 

The frequency- wavenumber spectra are shown in Figure 7(e)-(h). From the spectrum of the pristine plate, a clear 
A 0 mode can be seen. As in the v y results, additional wavenumber components (both positive and negative 
components representing forward and backward wave propagation) appear in the delamination cases. Figure 7(i)-(l) 
shows space-wavenumber spectra at 300 kHz obtained by short-space 2D FT. A window size of 8 mm, twice the 
wavelength of the A 0 mode in the full plate thickness, was used for the short-space 2D FT. As shown in Figure 7(1), 
the new wavenumber components only exist within the delamination region, which means this new wavenumber 
components are related to the waves propagating within the delamination region. The center of the new wavenumber 
component is around 5.6 rad/mm. However, this wavenumber does not match closely with the theoretical 
wavenumber (3.8 rad/mm) of the A 0 mode in the upmost two ply layers. The reason for this apparent mismatch is 
unclear. The A 0 mode is highly dispersive in the corresponding frequency-thickness range, and the bandwidth of the 
incident excitation could contribute to the mismatch. Another potential reason could be that the simulation spatial 
step size is not adequate for capturing the short wavelength behavior above the delamination (though recall that the 
step size is equal to 1 above the delamination). This unexpected mismatch will be further investigated in 

future work. 


CONCLUSIONS 

This paper presents our studies on guided wave propagation in CFRP plates using wavenumber-based data 
analysis methods. The frequency- wavenumber analysis enables determination of all existing wave modes while the 
space-wavenumber analysis reveals how the wavenumber component varies with respect to the propagation 
distance. 3D EFIT was used to simulate guided wave propagation in pristine and delaminated CFRP composites, 
generating time-space wavefields of the in-plane and out-of-plane velocity components. Waves are seen propagating 
within the delamination area at different velocities from the pristine case and in both forward and backward 
directions. 

Using the frequency-wavenumber and space-frequency-wavenumber analysis techniques, the wave interactions 
with delamination damage are more easily analyzed than when simply relying on time-space wavefields. Frequency- 
wavenumber results from both in-plane and out-of-plane motion studies show that there are new wavenumber 
components in the delaminated plate compared to those of the pristine plate. The space-wavenumber spectrum then 
reveals that the new wavenumber components appear within the delamination area. This finding is consistent with 
the standing waves observed in the wavefield images reported in the literature [13]. From the in-plane motion study, 
the new wavenumber component is consistent with the expected S 0 wavenumber for the upmost two ply layers. 
However, from the out-of-plane study, the new wavenumber component shows an apparent mismatch from the 
expected A 0 mode wavenumber in the upmost two ply layers. Further studies will be conducted to explore these 
results. 
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